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Abstract 

An efficient scheme for one-dimensional extensive air shower simulation and its implementation in the 
program CONEX are presented. Explicit Monte Carlo simulation of the high-energy part of hadronic and 
electromagnetic cascades in the atmosphere is combined with a numeric solution of cascade equations for 
smaller energy sub-showers to obtain accurate shower predictions. The developed scheme allows us to 
calculate not only observables related to the number of particles (shower size) but also ionization energy 
deposit profiles which are needed for the interpretation of data of experiments employing the fluorescence 
light technique. We discuss in detail the basic algorithms developed and illustrate the power of the method. 
It is shown that Monte Carlo, numerical, and hybrid air shower calculations give consistent results which 
agree very well with those obtained within the CORSIKA program. 

1 Introduction 

Monte Carlo (MC) simulation of extensive air showers (EAS) is the most common method to calculate de- 
tailed theoretical predictions needed for interpreting experimental data of air shower arrays or fluorescence 
light detectors. However, for primary particles of very high energy, straight-forward MC simulation is not a 
viable option because of the unreasonably large computing time required. The situation can be improved by 
applying some weighted sampling algorithms, like the so-called "thinning" method 1 1 1. i.e. treating explic- 
itly only a small portion of all shower particles and assigning each particle a corresponding weight factor. 
Although this approach allows the reduction of EAS calculation times to practically affordable values, it 
comes soon to its limits. The summation of particle contributions with very large weights creates signifi- 
cant artificial fluctuations for EAS observables of interest fTT'Tl. Imposing maximum weight limitations 
to ensure high simulation quaUty ||4J, on the other hand, prevents one from using less detailed sampling 
and correspondingly from further speeding-up the calculation process. A possible alternative procedure is 
to describe EAS development numerically, based on the solution of the corresponding cascade equations 
||5]|5|[71. Combining this with an explicit MC simulation of the most high-energy part of an air shower 
allows one to obtain accurate results both for average EAS characteristics and for their fluctuations (S). 

In this article we describe a new EAS simulation program of such a hybrid type, called CONEX. In 
CONEX the MC treatment of above-threshold particle cascading is realized in the standard way and does 
not differ significantly from the implementation in e.g. CORSIKA jO). On the other hand, the numerical 
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description of lower energy sub-cascades is based on the solution of hadronic cascade equations, using an 
updated algorithm of Ref. fTj, and a newly developed procedure for the solution of electro-magnetic (e/m) 
cascade equations. The corresponding algorithms are characterized by high efficiency and good accuracy 
even if a comparatively crude binning with respect to particle energy and depth position is used. Further- 
more, by accounting for neutrino production in addition to the typically considered particles, CONEX can 
also be used for the calculation of ionization energy deposit profiles. 

The numeric solution of cascade equations can be replaced, in principle, by a pre-tabulation of the 
characteristics of secondary sub-cascades, obtained via an iterative MC procedure IIOIII ll[T2l . An example 
of a combined approach, extended to a three-dimensional EAS simulation, is described in Ref. 1131 . where 
hadronic sub-cascades are treated using the method of Ref. Q and e/m sub-cascades are tabulated using 
the EGS4 MC program ifm . 

The approach presented in the current work does not require any pre-tabulation of particle cascades and 
is characterized by high efficiency and large flexibility. It can be applied to various initial conditions, i.e. 
a wide range of energies and angles of incidence of a primary particle, including the case of upward-going 
showers, as well as arbitrary parameterizations of the atmosphere of the Earth. These features make CONEX 
ideally suited for applications related to the event-by-event based analysis of EAS data, in particular, of 
fluorescence light based measurements. 

This is the first paper in a series in which we will investigate various features of EAS and their relation to 
the characteristics of hadronic multiparticle production. In this work, we shall present the hybrid simulation 
scheme in detail, leaving the study of shower predictions to a forthcoming article. 

The outline of the paper is as follows. Section|2]describes the calculation scheme and its basic proce- 
dures. In Sectionl^we show some examples of calculated EAS characteristics and investigate the accuracy 
and efficiency of the method comparing the hybrid approach with pure MC or numerical procedures. The 
reliability of the predictions is checked by a detailed comparison with CORSIKA results. Finally, a sum- 
mary is given in Sectionl^and both potential applications of the program and the prospects for its further 
development are discussed. 

2 Calculation scheme 
2.1 Physics overview 

The calculation scheme consists of two main stages: an explicit MC simulation of the cascade for particles 
with energies above some chosen threshold Ethr (being a free parameter of the scheme) and a solution of 
nuclear-electro-magnetic cascade equations for sub-cascades of smaller energies. Both MC and numerical 
parts are characterized by the same physics content, as described below. 

In the hadronic cascade one follows the propagation, interaction and decay (where applicable) of (anti-) 
nucleons, charged pions, charged and neutral kaons; all other types of hadrons produced in interactions and 
decays are assumed to decay immediately. Particle interactions in the MC part are treated within a chosen 
high energy hadronic interaction model (implemented are NEXUS 3.97 |Il21[l6|[ni, QGSJET 01 fT8l[T9l 
and II Il20ll21ll22l . and SIBYLL 2.1 1231 1241 l25l ). decays are simulated using the corresponding routines 
of the NEXUS model. The same models are used to pre-calculate secondary particle spectra for later use in 
the numerical treatment of hadronic cascade equations. Optionally, below some energy 100 GeV, 

GHEISHA 1261 is employed as low-energy hadronic interaction model. 

The MC treatment of the e/m cascade is realized by means of the EGS4 code 1 14|, supplemented by 
an account of the Landau-Pomeranchuk-Migdal (LPM) effect 1271 1281 l29l for ultra-high energy electrons 
(positrons) and photons. The simulation of photonuclear interactions and muon pair production was added 
to the EGS package closely following the CORSIKA implementation f9l. The system of coupled e/m cascade 
equations is based on the same interaction processes as implemented in the MC, using Bethe-Heitler cross 
sections for the bremsstrahlung and pair production with energy-dependent correction factors according 
to Storm and Israel 1301 . the Klein-Nishina formula for the Compton process, and accounting for Moeller 
and Bhabha processes as well as for positron-electron annihilation (see, for example 1141 1311 ). Both LPM 
suppression and photo-effect are neglected in the e/m cascade Eqs., as the latter are employed in the energy 
ranges where these processes are not important. Ionization losses of electrons and positrons are described 
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by the Bethe-Bloch formula with corrections to account for the density effect 1321 . 

High energy interactions of muons (bremsstrahlung, pair production and muon-nuclear interaction 1331 
I34II35I 36 371) are taken into account in the MC part but are neglected in the cascade Eqs. 

In general, an individual shower is simulated as follows. One starts with the primary particle of given 
energy, direction and initial position in the atmosphere (by default, at 100 km above sea level, if no special 
geometry is required, e.g. up-going showers). The initial particle direction thus defines the position of 
the "true" shower axis, in the following referred to as the "shower trajectory". For a hadron as primary 
particle, one simulates the hadronic cascade explicitly, recording all secondary particles at a number of 
pre-chosen depth levels and energy intervals, until all produced secondaries have an energy lower than the 
threshold Ethr- The levels are defined with respect to the projected depth X, i.e. the slant depth for the 
particle position projected to the initial shower axis (shower trajectory), as described in more detail in the 
Appendix. All sub-threshold hadrons/muons and e/m particles are filled into energy-depth tables that form 
the "source terms" for the cascade equations. In parallel, the above-threshold e/m particles are transferred 
to EGS4 for simulating the e/m particle cascade in a similar way, with all sub-threshold e/m particles being 
added to the e/m source terms. 

In the next step the hadronic cascade at energies below -Ethr is calculated numerically for the first depth 
level using the corresponding cascade equations and initial conditions specified by the source terms. As 
the result, one obtains discretized energy spectra of hadrons of different types at the next depth level. All 
sub-threshold e/m particles produced at this stage are added to the e/m source term. Then sub-threshold e/m 
cascades are calculated by solving the corresponding e/m cascade equations for the given initial conditions. 
Hadrons due to photonuclear interaction and pair-produced muons that are generated in the numerical 
solution of the e/m cascade Eqs. are added to the hadronic source term of the next slant depth level. This 
procedure is repeated for the following depth levels, each time using the hadronic and e/m source terms of 
the previous level. 

Ultra-high energy e/m particles can undergo geomagnetic pair production and bremsstrahlung well 
above the atmosphere of the Earth I38ll39ll40l . Therefore, in case of the primary particle being a photon or 
an electron, the simulation process starts with the calculation of possible interactions with the geomagnetic 
field using the PRESHOWER code 1411 and the above described procedure is applied to the secondary 
particles. 



2.2 Hadronic cascade equations 

The backbone of a hadron-initiated extensive air shower is the hadronic cascade which develops via particle 
propagation, decay, and interaction with air nuclei of both the initial particle and of produced secondary 
hadrons. The corresponding integro-differential equations are given by |7| (see also (01) 

dh^iE,X)\, h^iE,X)\, H^^E,X)\,]%- + 4-{priE)KiE^X)\,) 



dX XaiE) ' Ta{E)c dE 

+ V / dE' hd{E\X)\j, 

., JE 



d 

-fhad 



Wd^a{E',E) |#|^ • 

—^^^(^ + D,^a{E.E)-^^ 



+ (1) 

where ha{E, X)\rp are the differential energy spectra of hadrons of type a with energy E at depth position 
X along a given straight line trajectory T (in the following the T-symbol will be omitted), /3^°"(i?) = 
~dEa/dX is the ionization energy loss of particle a per depth unit. A muon is treated like a hadron, but 
without interaction term. 

The first term in the rh.s. of Eq. Q represents the decrease of hadron number due to interactions with 
air nuclei 

dhg _ ha 
dX~ A,' 

with the corresponding mean free path Xa = n^aii /f|^or" ' where mail is the average mass of air molecules 
and c°i^f is the hadron a - air nucleus inelastic cross section. 
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The second term describes particle decay, with the decay rate on a path dL being 

dha = -ha—, (3) 

Ta C 

where Ta is the life time of hadron a in the lab. system, related to the proper life time ri"^ by Tq = 
ri"' E/rria, with rria being the hadron mass and c the velocity of light. From the definition of slant depth 
i28t follows 



dL 



dX 



(4) 



Pair(X) 

The third term in Eq. Q takes into account particle ionization energy losses and the integral term in 
Eq. Q represents the production of particles of type a in interactions and decays of higher energy parents 
of type d, with Wd^a{E\ E), D^^aiE' , E) being the corresponding inclusive spectra of secondaries. 

Finally, the so-called source term 5)j""^(£',X) defines the initial conditions and is determined during 
the MC simulation of above-threshold particle cascading. It consists of contributions of all sub-threshold 
hadrons produced at that stage 



= 5fc^had(£;^^) ^ ^ 5-^S{^E~E,)5{X~Xi), (5) 



i=l 



with di, Ei, Xi being type, energy, and depth position of the source particles. 

The numerical method of solving the hadronic cascade equations is similar to the approach of (TJ and 
is summarized in Appendix l5.2l 

2.3 Electro-magnetic cascade equations 

The e/m cascade development can be described by the following system of integro-differential equations 
(see, for example, 1311 ') 



^^^-^^'^^ -aAE)lAE,X) + ^{l3i°2\E)k-iE,X)) 



dX c y , c y . , ■ 



f '^^^dE' [l,-{E',X)W,-^,-{E',E) + l,+{E',X) 

J E 



X W,+_,AE\ E) + l^E', X) W^^AE', E)] + S^L'^E, X) (6) 

= -aME)lAE,X) + ^{piriE)l,,iE,X)) 

+ / dE' [lAE',X)W,+ ^e+{E',E) 
Je 

+ 1,{E\ X) W^^AE'. E)] + S^jAiE, X) (7) 

^^^^^ = -a,{E)l,{E,X) + j^J'^^dE' 

X [1^{E', X)W^^iE', E) + /,-(£;', X) W^-_^{E', E) 

+ lAE', X) W,+_iE', E)] + S;^^E, X), (8) 

where la{E, X) (a = e~ , e+ , 7) are energy spectra of electrons, positrons, and photons at depth' X, <Ja{E) 
are interaction cross sections (in units area/mass, see Sec. 12. 2> 



~ "'(brcmsstrahlung) ~r '''(Mocllcr) 

'-^(bremsstrahlung) ~^ ^(Bhabha) ~t~ ^^(annihilation) 



"'7 — "'(pair production) + O'(Compton) + C(photonuclear) + C'(,nuon pair): (9) 



' In the absence of particle decays there is no dependence on a particular shower trajectory, apart from the density effect correction. 
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and Wd^a{E' , E) are corresponding differential energy spectra of secondary particles 



W^^_^AE', E) = Ty]'"^|-(S', E) + W^^-°^l"(E', E) 

W,+ _,,+{E\ E) = iy^^™;^(S', E) + WfJ^$^+%E', E) 

W^-.^{E\ E) = W^^'^'^'°\E\E) 

W,-^^{E', E) = W^l^Zl-iE', E' - E) 

W,+ _,^{E', E) = W^l^^l-{E\ E' -E) + W^^'f^iE', E) 

We+-,e-iE', E) = Wf:^$^+\E', E' - E) 

W^^^E', E) = W'^Z\-{E',E) + W^^^^°\E', E' - E) 

W^_,AE', E) = W'^ZViE'^ E). (10) 

Here cr(yiac\\cr) ™d o'(Bhabha) Correspond to the process of (5-electron knock out above some energy thresh- 

MocUcr/Bhabha/ ,-,\ / 77-^/ TirMocllcr/Bhabha/ 7-, 7-n/\ / 1 1 \ 

<^e-/+ {E)^ dEW^_^^J^_^^ (.E,E), (11) 

whereas the contribution of those processes below i?"/™ is treated as continuous energy losses and consti- 
tutes a part of PfoJ.E) 

dE,± 



mE) 



dX 



ionization+(5(<£J^." 



(12) 



The bremsstrahlung cross section diverges due to the characteristic infra-red singular behavior l/£" of 
the secondary photon spectrum W^1°^^{E, E') = E — E') and normaly requires to introduce 

some low energy cutoff. The sub-cutoff photon emission could be treated as continuous "radiation" energy 
losses as is done in EGS4 il4l . However, this is not necessary in our calculation scheme as is shown in 
ADPendix l5.3l where the numerical solution is presented. 



2.4 Photonuclear effect and muon pair production 

To take into account hadron and muon production by photons, two additional terms are introduced that 
couple e/m and hadronic cascade equations via source terms. Photoproduction of hadrons, i.e. photonu- 
clear interaction, is implemented using the cross section of Ref. 1421 . Particle production distributions are 
approximated by those of Tr'^-air interactions. Then the corresponding source term can be written as 

^.cnwhad^^;^^) ^ / dE' l^{E',X)W^^^^a{E',E)af°'°'''"'{E'), (13) 

JE 

where W.,,o^a{E' tE) is defined in analogy to Eq. (I10> . With a small cross section, a photon can also 
produce a yU^/i^ pair This gives another contribution to the hadronic source term 

S-J-^t^{E,X)= / dE' l^{E',X)W^^^{E\E)a"^'\E'). (14) 

JE 

Finally the source term defined in Eq. (|5|i becomes 

S'^'^'^iE, X) = Sf^-^^'-^^iE, X) + S^J^-^^'-^^iE, X) + S'™^^'^(£:, X) . (15) 



3 Applications 

In the following we demonstrate the reliability of the CONEX code by comparing predictions for shower 
observables calculated with cascade Eqs. and the full hybrid scheme to that of MC simulations. As CONEX 
can also run in pure MC mode, both CONEX and CORSIKA are used for calculating the MC predictions. 
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Figure 1: Average hadronic shower size profiles (left panel) and energy spectra at X = 500 g/cm^ (right 
panel) of nucleons and charged pions for proton-initiated vertical {0 = 0°) showers of 10^^ eV. Compared 
are the results of solving numerically the system of cascade equations (CE) with different discretization 
bin sizes in energy and depth. 



If not specified otherwise, the calculations were performed using the QGSJET model for hadronic in- 
teractions at energies E > ijjj,^ = 80 GeV and GHEISHA model for E < E^^. Our default choice for 
the cutoff between the MC and the numerical parts is i?thr = lO^^-Eo, Eq being the energy of the primary 
particle. The default energy grid for solving the cascade Eqs. is 30 bins per energy decade (ds ~ 30) for 
hadrons and e/m particles and the slant depth binning has a 5 g/cm^ elementary step (AX). When apply- 
ing the hybrid scheme, high energy particles are treated in MC. As a consequence the energy transfer from 
hadronic to e/m particles is more precise, allowing us to use larger bins, i.e. 20 bins per energy decade and 
a 10 g/cm^ slant depth step size. 

3.1 Hadronic shower component 

In Fig.^we investigate the stability of our scheme and compare both longitudinal profiles of nucleons and 
charged pions and their energy spectra at 500 g/crn^ for different choices of energy and depth discretization. 
Results only change significantly for very large discretization intervals. 

In Fig.|2lwe plot similar characteristics of charged pions and muons for 10^^ eV proton-initiated show- 
ers simulated with QGSJET 01 at high energy and GHEISHA at low energy. The results are compared to 
CORSIKA predictions. The agreement between the results from the different CONEX calculation methods 
as well as CORSIKA simulations is very good. 
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Figure 2: Average longitudinal shower size profiles (left panel) and energy spectra (right panel) of charged 
pions and muons with energies above 1 GeV. The calculations were done for proton-initiated vertical show- 
ers of 10^® eV. Compared are the predictions obtained with CONEX applying the hybrid (dashed line), pure 
MC (points) and numerical calculation (full line) schemes. In addition CORSIKA predictions are shown 
as symbols (stars). 

3.2 Electromagnetic shower component 

The longitudinal profiles of electrons, positrons, and photons for a 10^'* eV vertical photon-initiated shower 
are shown in Fig.|3(left panel). The shower size profiles are given for the cutoff energies i?^;™ =1 MeV, 
1000 MeV using again the hybrid, pure MC, and cascade equation approaches. 
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Figure 3: Average longitudinal shower size profiles of charged particles for £','^/™ =1 MeV (top left panel) 
and E'^^^^ =1000 MeV (bottom left panel). Particle energy spectra of photons, electrons, and positrons 
are given in the right panel. All curves are calculated for photon-initiated e/m showers using both hybrid 
(dashed line), pure MC (points) and cascade equation (solid line) approaches. 

While for a large cutoff energy, for example 1000 MeV, the agreement between the different meth- 
ods is good we notice systematically larger particle numbers in the hybrid and numerical calculations for 
^min ^ 1 MeV. The corresponding difference is clearly visible in the particle energy spectra and is related 
to spatial effects in the shower development. Low energy electrons (positrons) undergo significant angular 
deflections due mainly to multiple Coulomb scattering. In turn low energy bremsstrahlung photons pro- 
duced by such deflected particles also have significant directional deviations from the initial shower axis. 
If only the track projected to the shower axis is considered, this leads to an apparently faster absorption 
of low energy particles (higher interaction rate and ionization energy loss) compared to that expected for 
particles traveling along the shower axis only (see also discussion in 1431 ). Although a full account of this 
effect requires a three-dimensional treatment of the particle cascade at MeV energies a reasonable improve- 
ment can be achieved by introducing an "average angular deflection". As the effect is only important for 
low energy leptons which anyway lose their energy quite fast we may estimate the corresponding average 
scattering angle of an electron (positron) as 

{e') ^ ^m. (16) 

where Es — 21 MeV and L{E) is the average travel distance of an electron of energy E in units of radiation 
length. With L{E) ~ E/E„it and E^.a ^ 81 MeV, we have 

{e^) ~ . (17) 

^crit 

For numerical applications, expression Ml\ has to be modified to satisfy the boundary condition < 7r/2. 
In the following we chose the ansatz 



1 — cxp 



E^ 
E 



(18) 
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which reduces to the functional form of for large E and approaches (7r/2)^ in the small E limit. The 
same kind of formula is used for photons but with a different parameter E!^^, as photons themselves do 
not undergo multiple scattering. Good agreement between fully three-dimensional MC simulations and 
cascade equation calculations is obtained for E°± = 9.5 • lO"** GeV and E°^ = 5 • 10"^ GeV. This is 
shown in Fig.|4|where the shower size profiles and energy spectra of the two approaches are compared for 
e/m showers of 10^^ eV and 10^^ eV. 




depth (g/cm-) energy (GeV) 

Figure 4: Left panel: Average longitudinal shower size profiles of charged particles for 10^^ eV (top) and 
10^^ eV (bottom) photon-initiated showers and S^/j™ = 1 MeV. Right panel: Particle energy spectra at 
X ~ 700 g/cm^. Shown are the results of MC simulations with CONEX and CORSIKA and of the hybrid 
and cascade equation schemes. Here an "average angular deflection" is used in the numerical and hybrid 
schemes. 
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Figure 5: Left panel: Average longitudinal profiles of charged particles and photons of energies above 
1 MeV for proton-initiated vertical showers of Eq = 10^^ eV. Right panel: Particle energy spectra of 
photons, electrons, and positrons for the atm. depths X = 700 and 1000 g/cnP. Shown are the results from 
the hybrid calculation (dashed line), pure MC simulation (points), and numerical cascade Eqs. solution 
(full line). In addition CORSIKA predictions are given by stars. 

After having shown that both hadronic and e/m showers can be well described by the presented cas- 
cade equations and the corresponding hybrid simulation scheme, we test the coupling between the e/m 
and hadronic shower components. In Fig.|5]we compare both longitudinal profiles and energy spectra of 
e/m particles for vertical proton-induced showers of 10^^ eV. Very good agreement is found between the 
results of the different calculation methods and also with the CORSIKA predictions. The simulations were 
performed using QGSJET 01 as high-energy interaction model and GHEISHA for hadronic interactions at 
low energy. 

3.3 Energy deposit 

The longitudinal profile of simulated EAS depends on the e/m low-energy cut-off E^^(-^ used in the calcu- 
lation (in this work 1 MeV). Lowering this threshold to, for example, 50 keV would increase the number 
of charged particles at the shower maximum by a few percent ||44|. In the case of the photon number, 
the dependence of the simulation results on the low-energy cut-off is much stronger. In fact, the number 
of photons diverges if the simulation cut-off is set to 0. Problems of this kind can be avoided if, instead 
of the secondary particle profile, the ionization energy deposit profile of showers is considered. Another 
advantage of the energy deposit profile is its direct relation to the light curve measured in air fluorescence 
experiments. Current measurements support the theoretical expectation that the fluorescence yield is pro- 
portional to the ionization energy deposit 1451 . 

Therefore, the hybrid simulation scheme implemented in CONEX has been extended to also allow the 
calculation of energy deposit profiles as described below. 

The calculation of the energy deposit profile in the MC part of CONEX is very similar to that in COR- 
SIKA, see 1461 . The deposited ionization energy is counted for each particle and traversed slant depth bin. 
If a particle reaches the cut-off energy of the calculation, depending on the particle type, either all its en- 
ergy or a part of it is assumed to be deposited locally. For example, a positron is expected to annihilate in 
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Table 1: Energy deposited locally for particles with £'ki„ < i?min (nic denotes the electron mass and niM 
is the nucleon mass). 



particle 


7 


e 


e+ 
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Etot 
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Etot — "ijv 


Etot + riiN 


Etot - A ■ niN 


Etot + A ■ rriN 



the end and, therefore, all its energy plus electron mass are assumed to be deposited. The different energy 
contributions taken to be deposited locally are given in Table[T] 

Solving cascade equations numerically one can calculate the energy deposit of particles with E > E'min 
explicitly. In addition the number of newly created particles at each step in atmospheric depth is given by 
the source functions. However, the energy carried by particles falling below the cut-off energy threshold 
and that of neutrinos is not directly calculated. 

For the e/m cascade equations, the total energy deposit per slant depth bin can be estimated from energy 
conservation. The deposited energy follows from the difference between the total kinetic energy at slant 
depths Xjn and including a correction for positrons similar to the MC treatment 



Eii;ix„ 



(19) 



where me is electron mass. 

A similar energy conservation-based method can be applied to the hadronic shower component, how- 
ever, the deposited energy and the energy going into neutrino production have to be distinguished. Since 
we know the ionization energy deposited by hadrons and can calculate the number of neutrinos produced 
within a given slant depth bin, we can obtain the energy of all the particles falling below the low-energy 
threshold from 

Ecut{Xm) ~ E]:,al{Xrii) — Eion{Xrii) — E^{Xm), 



where 



E[on (^m) 



£ iEf^^Y.''aiXm)+J2maK{X„,)\ 

— imin \ CL / 

E ( E K{x^-i) + E "^'^ K{x„,-i) 
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X exp 
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had\ 
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/Olon f 77'had^ 
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^had ^had 

rUa \L{Xm) - L{X,n^i) 



E^{Xjn) 



^^(O)^had 



d — 



d j=i+l 



P.UX') 



(20) 



(21) 



(22) 
(23) 



To account for the fact that only a part of the energy carried by hadrons falling below the low-energy 
threshold is deposited as ionization energy, a factor /had//^ = 0-45 is introduced f47J . Thus, the energy 
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deposit for the hadronic part of the system of cascade Eqs. is given by 



1 



(24) 



Summing the different contributions, the longitudinal energy deposit profile can be calculated. A com- 
parison of the results of different calculation methods within CONEX with CORSIKA predictions is shown 
in Fig.|6l(top panel). Good agreement is found. 

An interesting consistency check of the energy deposit calculation scheme is the investigation of the 
dependence of corresponding results on the low-energy threshold used in the e/m cascade Eqs. In Fig. |6l 
(bottom panel), the mean longitudinal energy deposit profiles of iron-initiated showers of 6* = 60° and 
E = 10^^ eV are shown for different low-energy thresholds. The profiles are independent of this threshold 
over a wide energy range - the approximation of local energy deposit breaks down only for a threshold 
energy of 100 MeV and higher. 
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Figure 6: Upper panel: Average longitudinal energy deposit profile for vertical proton-initiated showers of 
10^* eV. The results obtained with the e/m energy cutoff E';^.'^ = 1 MeV, using hybrid (dashed Une), MC 
(points), or numerical approaches (full line), are compared to CORSIKA predictions (stars). Lower panel: 
Energy deposit profiles for 10^*^ eV inclined [6 = 60°) iron-initiated showers for the e/m energy cutoff 
=1 (full Hne), 5 (dashed line), 10 (dotted fine) and 100 MeV (dashed-dotted line). 

3.4 Shower fluctuations 

So far we have only studied mean shower observables, i.e. ones averaged over many showers. As a proper 
description of shower fluctuations is of central importance for the analysis of experimental data, we will 
also discuss the treatment of shower fluctuations. 
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Running CONEX in hybrid mode allows us to benefit from the fast numerical solution of cascade Eqs. 
and, at the same time, to obtain a good description of shower-to-shower fluctuations. Here, the key param- 
eter of the method is the energy threshold that separates the explicit MC simulation from the application 
of cascade Eqs. By default, this energy threshold is set to Ethi = 0.01 • Eq for all particles in CONEX. In 
principle, it can be chosen differently for e/m and hadronic particles to further reduce the simulation time 
needed for high-energy showers. 




Figure 7: Fluctuations of the shower maximum depth XmaxSround the mean shower maximum depth 
(X,„ax) for a primary energy of lO^^eV for proton and iron-initiated showers simulated with CONEX (lines) 
and compared with CORSIKA results (points). 

In Fig. Q the distribution of the depth of shower maximum is shown for proton- and iron-initiated 
showers. Within the statistical uncertainties, the CONEX results agree very well with that obtained with 
CORSIKA MC simulations. Not only the fluctuations but also the mean depth of shower maximum obtained 
with the two codes are the same (not shown). 

The dependence of the simulated fluctuations on the energy threshold is shown in Fig. |8] Even a 
threshold as large as Ethi = 0.99 • Eq is sufficient to reproduce almost the full X^ax distribution. This 
comparison demonstrates that the fluctuations of particle production in the first interaction of such high- 
energy showers determine almost the entire shower profile. 

4 Summary 

We have developed a fast and efficient one-dimensional hybrid simulation scheme for ultra-high energy air 
showers. It combines explicit MC simulation of high-energy particle interaction, propagation and decay 
with the numerical solution of a system of cascade equations for calculating the low-energy part of the 
particle cascade. 

The presented hybrid simulation scheme is implemented in the code CONEX^. Several high- and low- 
energy hadronic interaction models are available within CONEX to study theoretical predictions and the 
model-dependence of data analyses. 

All relevant interaction and decay processes are considered in both the MC and the cascade equation 
parts of CONEX. These processes also include muon pair production and photonuclear interactions of 
muons. At ultra-high energy, the LPM effect and possible e/m preshowering in the geo-magnetic field are 
simulated. 

^The CONEX code is available upon request from Tanguy.Pierog@ik.fzk.de 
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Figure 8: Fluctuations of the shower maximum depth X^ax for lO^^eV proton-initiated showers simulated 
with CONEX for three different energy thresholds Ethr/En = 0.01, 0.1, 0.99. 

The hybrid simulation scheme has been extended to include the simultaneous calculation of both shower 
size profiles of various particles and the generation of ionization energy deposit profiles. The latter are 
independent of the low-energy cut-off that has to be applied in all shower simulations. Knowing both 
the shower size profile (with an arbitrary low-energy cut-off) and the energy deposit profile allows us 
to simulate directly the fluorescence and Cherenkov light signal of air showers. Together with the fully 
three-dimensional implementation of the shower axis geometry, this makes CONEX ideally suited for event 
simulation and data analysis of fluorescence light experiments such as HiRes f481. Auger 1491 . TA 1501 . 
and EUSO 1.51 J . 

In developing CONEX, particular emphasis is put on the accuracy and reliability of the shower simula- 
tion to make the code directly applicable to data analysis of air shower experiments. Extensive comparisons 
with CORSIKA simulations show that all shower distributions agree very well. Both mean shower profiles 
and energy distributions as well as their fluctuations were compared, only a small fraction of which could 
be shown in this paper. 

In a forthcoming work we will study the influence of different hadronic interaction models on air 
shower predictions. In particular we will investigate the total calorimetric energy deposited by a shower in 
air First results of this work have been presented in 1521 . 

5 Appendix 
5.1 Geometry 

Using a one-dimensional treatment of EAS development, a given shower trajectory may be characterized 
by a single parameter - its distance to the center of the Earth O: i?^ = - see Fig.|9l where Pf 

is the lowest trajectory point. In case the observer is positioned at some height /lobs above sea level and 
^Earth + ^obs > R± (^^Earth being the Earth radius) the observed shower inclination is 

9 = arcsin — — . (25) 

^Earth + Alobs 

The position of a particle moving along a given trajectory may be then characterized by its local az- 
imuthal angle 6p with respect to the vertical direction, or alternatively, by its height Hp, or by the distance 
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Figure 9: Geometry for defining the shower trajectory accounting for the curvature of the Earth. 



Lp ~ IPiPl to the lowest trajectory point 



- 

Lp — 



smf 



cos 6*1^ 



(26) 
(27) 



For the solution of cascade equations, it is more convenient to use instead the slant grammage Xp, i.e. 
the integral over the atmospheric density pnu {h) along the trajectory T from the point P to infinity 



dl' Pair {H{1')) 



(28) 



and characterize an arbitrary particle position by two variables, and X. 

5.2 Numerical treatment of hadronic cascade equations 

In order to reduce the problem to the solution of a system of ordinary differential equations, we perform 
a discretization of particle energy spectra. Introducing an energy binning as Ef^'^ = -E-mtn C'^^ C ~ 
lO^/d-E (ijhad ^ I ds = W ^ 30) and replacing the smooth particle spectra ha{E, X) by discrete 

*max) we may get instead of ([T) 

Pa ) 



contributions h^J^X) of representative particles of energies E^ {i = 1 



dKjX) 
dX 



-Kix) 



dL 



dX 



d j=i 



+K^\x) 



d^a phad 



TTihad TTihad 



(29) 



where we used Q, replaced the integral over parent particle energies / dE' by the discrete sum 
introduced the discretized source term S]^f'^{X) and the discrete interaction (decay) spectra W^i^^ (^d^o) 
via 



max[£;i';'^,£;;^';^] 



dE' S^'"^{E',X) KiE'/Ef""^) 
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min[E]'f^,Ethr] 



had 



_ghad 



dE' S^^HE\X) [1 - K{E'/E]^f,)] 



dE' Wd^a{E' , Ef""^) KiE'/Ef""^) 



dE' Wd^aiE'^E'l^") [1 - KiE'/E^f)] 



(30) 



(31) 



and the discrete energy loss term via : 



/Dion / 77'had\ 



KiX) 



^had _ ^had 



Here the condition of both energy and particle number conservation gives 

Ce-1 



K{e) = 



C- 1 



The solution of the homogeneous part of Eq. M9\ is 



h\{X) = K{X^) cxp 



1 - VF^^, 



Pa y^i 



^had _ ^had 



{X - Xo) 



ma\L{X) - L{Xo)\ 



(32) 



(33) 



Correspondingly the solution of the full equation can be obtained in an iterative way 



K{x) 



KiXo) exp 



XaiEf^'^) £:had_£:had 



{X - Xo) - 



ma\L{X)^L{Xo)\ 

„(0) phad 



oion / j^hsbd \ 



^had ^liad 



X exp 



1 - 



had^ 



\d{E^ 

^had ^had 



had\ 



mdi{cn ) 



d^a pihad 



ma\L{X) ~ L{X')\ 



(34) 



Discretizing depth positions as X^ = mAX, m = 1, ...,m,„ax, -''^m.nax corresponding to the ob- 
servation level for the given shower trajectory, the formulas ( I33l34t can be used to calculate spectra of 
different hadrons at all depths Xm, starting from the initial condition /i^(Xm) for m = 1, calculating first 
using Eq. 03, and h];;^-'^-^{X„i+i), ... , hl{Xm+i), using Eq. (|34), etc. 

for m = 2, TTiniax- The depth integral in Eq. ( I34> is taken using the Simpson formula; particle spec- 
tra values /Iq(X„j_|_i/2) mids of the bins = (-^^m + Xjn+i)/2 are obtained via a logarithmic 
interpolation between previously calculated hl^{X„i), /i^(X,„+i). 

As for neutral pions, we assume them to decay at the place, calculate the number of tt^s produced at 
depth X„i as 



X„ 



dx'Y: E ^'.(^') 



d j=i+l 

ji rnd/icTf^) 



D 



d—*Tr° T^had 



Pair(^') 



X„ 



dx' sy^{x') 



(35) 



and add photons resulting from 7r° decay to the e/m source function. 
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5.3 Numerical treatment of e/m cascade equations 

To reduce the problem to the solution of a system of ordinary differential equations we again, like in case 
of hadron cascade, discretize particle energy spectra la{E, X) using an energy grid £'°''™ = ^niin C**^ 
(E'^-^ = 0.1-^1 MeV), replace the integral over parent particle energy E' in ( I6I8> by a sum over discrete 
energies, and discretize the source term Sa^"\E, X) and the differential energy spectra Wd^a{E' , E) 
according to (QUI?!! (with Ef""^, E^^^"^, S^^'^E' , X) being replaced by E"/"", E"^-^, Sl'^^iE' , X)). Then 
Eqs. ( I6I8> are transformed to 

^ - E E I'M) + Sli-iX), (36) 

d j=i 

where we replaced the continuous energy loss term {^f3^°±{E)l^± {E, X)) by 



-Pe±\E, ) + li^± i-t'i+i ) —j^ — —j-^ 



J^iJ/ 111 _ ^U/ 111 ^ - . , ^ . ,u ^ 

and included the two terms, together with the interaction cross sections, into the discrete particle production 
spectra W^^_^^, defined as follows 

m^e^ = - aAEl'-) - ^'1 (37) 

and W'll_^^ = VF^^jj for all other combinations of i, j, a, d. 

It is worth verifying that all elements of the matrixes Wj'^^ are free of singularities. Indeed, singular 

terms appear in a^'^^'^'XE'l'''^) = Sf''"'dE' W'^^'^l^iEl' E') and in W%_^^± (defined by Eq. ^) 
due to the characteristic \/{E — E') dependence of W^l°™l_[E, E'). Nevertheless, in W^±^^± defined 
by Eq. ( I37> the corresponding singularities are canceled against each other 



W^e±-le± = W"^^"""""^' -cr°i™''(i;^' ) (40) 

c/„. CE'/E:/'^^-! 



f ' dE' W^I''^!-{Ef"\E') 



C-l 



rf^' iybfcms_(£,c/m^^,-| (41) 







C 



e:/'-(c-i)Je:^ 



rEl'"' 

/ dE' (£;°/'" - E') W];i''^l-{Ef'\ E') 



where the last two integrals are finite. Similarly one can check the finiteness of W^±^ ± . 
To solve the system ( I36> . we first find the solution of the homogeneous equation system 



^^=E^r-.a^^W, (43) 
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which is given as 

3 

i:{X) = Y,D:j{Xo)e^'^''-'''>^ (44) 



1=1 



Here A} and DI^j{Xq) are correspondingly the eigenvalues and eigenvectors of the system of linear alge- 
braic equations obtained by inserting ( I44> into ( I43> : 

d 

with the normalization fixed by the initial conditions 2S.X = Xq: Y^'j=i ^ai i-^o) ~ 'a(^o)- 
Then the solution of the full equation system i36\ may be given in a recursive form 



1=1 



gAi (x-x„) ^ J2 F}^{Xo) fdX' 

d -^^0 



(45) 



j=i+l g 

where the matrix Fj^ is inverse with respect to D^^j 

J2DUXo)FUXo)=Si. (46) 

Equations ( I44l45t can be used to calculate the discrete spectra of e/m particles (ATm) at all depths 
Xm in the same way as in case of the hadronic cascade. 
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